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We discuss the relevance of directional detection experiments in the post-discovery 
era and propose a method to extract the local dark matter phase space distribution 
from directional data. The first feature of this method is a parameterization of 
the dark matter distribution function in terms of integrals of motion, which can 
be analytically extended to infer properties of the global distribution if certain 
equilibrium conditions hold. The second feature of our method is a decomposition 
of the distribution function in moments of a model independent basis, with minimal 
reliance on the ansatz for its functional form. We illustrate our method using the 
Via Lactea II N-body simulation as well as an analytical model for the dark matter 
halo. We conclude that 0(1000) events are necessary to measure deviations from 
the Standard Halo Model and constrain or measure the presence of anisotropics. 
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I. INTRODUCTION 

Dark matter comprises 80% of the mass of the Universe and its presence has been inferred 
gravitationally in numerous different physical settings including galactic rotation curves, 
precision studies of the cosmic microwave background radiation, and through gravitational 
lensing of galaxies and clusters of galaxies. However, dark matter's non-gravitational 
interactions remain unknown and there is an active search underway to determine dark 
matter's identity. 
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Many well motivated proposals to extend the Standard Model accommodate naturally a 
dark matter candidate that is stable, electrically neutral and has the correct relic abundance 
PP . These models generically predict that dark matter should interact with Standard Model 
particles with electroweak-sized cross sections. Currently, a large experimental effort is 
underway searching for interactions of dark matter with nuclei in underground experiments. 

Typically, direct detection experiments attempt to infer that a dark matter particle 
scattered off a nucleus by detecting the nuclear recoil and measuring the energy deposited 
in the detector. While most presently running experiments were designed to only measure 
the energy of the nuclear recoil, more can be learned about dark matter if the direction 
of nuclear recoils is measured as well [2j. The anisotropy of dark matter-induced nuclear 
recoils partially arises from the motion of the Earth through galactic dark matter halo, but 
can also arise from intrinsic anisotropics in the dark matter halo [3j. Therefore, measuring 
the direction of the nuclear recoils is not only a powerful handle to discriminate signal 
from background but also provides a whole class of new measurements of the dark matter's 
velocity distribution in the solar neighborhood [3HH]. 

Over the past decade, the first experiments capable of detecting the directionality of 
nuclear recoils have been built and commissioned [X0] - fl9] . More recently, the first prototypes 
have started to release results that are approaching 0.1 kg-day exposures [201423] . 

The majority of the upcoming directional detection experiments use time projection 
chambers (TPCs) to resolve the direction of the nuclear recoils, by means of drifting the 
ionized track originated from the nuclear recoils through an electric field and projecting it 
over a charge collection grid. Two of the dimensions are then inferred from this projection, 
while the third dimension is obtained by the drift time of the particles in the ionized track. 
In order to resolve the direction of the track to good accuracy, these TPCs have to operate 
at low pressure and density, and therefore are limited in target mass, as well as in volume, 
since the drift length is limited by diffusion [T3"l \W\ . The limitation in exposure of directional 
detection experiments makes them currently not competitive with current experiments such 
as XENON100 [21] and CDMS [25], amongst many others. 

Most likely, directional detection experiments will not be able to scale their exposures 
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sufficiently fast to be the first to discover dark matter. However, they will provide unique 
measurements of dark matter properties, once dark matter is discovered through direct 
detection. The unique capabilities of directional detection experiments have been under- 
emphasized. For instance, directional detection experiments will explore the kinematics of 
the interaction between dark matter and Standard Model nuclei more directly than direct 
detection experiments [26143T] . 

This article will explore how directional detection experiments can measure astrophysical 
properties of dark matter, such as its phase space distribution in the solar neighborhood. 
Previous studies in this domain have been carried out in [32-36J assuming specific 
parameterizations for the functional form of the velocity distribution. One important goal 
of this article is to investigate the dependence of the dark matter phase space distribution 
in different coordinates, such as integrals of motion. If the dark matter distribution in the 
solar neighborhood is dominated by the equilibrium component, measuring the dependence 
of the local distribution on integrals of motion will allow us to infer properties of the global 
distribution. Another goal of this article is to remain agnostic about the functional form 
of the local distribution and explore generic parameterizations in terms of a basis of special 
functions. 

More properties of the entire galactic halo can be inferred by combining the measurement 
of the phase space distribution function with N-body simulations, ultimately giving clues to 
the original formation of the Milky Way galaxy. 

The organization of this article is as follows. In Sec. [IT] we overview general aspects of 



directional detection in the context of post-discovery. In Sec. |III| we discuss Jeans theorem 
and how it motivates integrals of motion as powerful variables in terms of which the dark 



matter distribution function can be expressed. In Sec. |IV| we present a novel method to 
parameterize the functional form of the distribution function - instead of relying on ansatzes 
of its functional form, we perform a special function decomposition of the dark matter 
phase space distribution. The coefficients of this decomposition will be the measurable 
quantities in a directional experiment. We perform fits of these expansion coefficients in 
Sec. IVlby simulating data from a mock experiment, and discuss the level of precision in the 
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measurement the velocity distribution function that can be achieved with plausible numbers 



of events at upcoming directional detection experiments. In Sec. VI we test our method using 
the Via Lactea II N-body simulation as a template of the dark-matter distribution in the 
solar neighborhood, by scattering the particles of the simulation at a detector placed at 8kpc 
from the center. We then use this mock signal to recover the input distribution, providing 



a closure/bias test for the procedure. Finally, in Sec. VII we discuss the implications for 
upcoming experiments. 

II. GENERAL ASPECTS OF DIRECTIONALITY 

In this section we discuss how the direction of nuclear recoils captures the underlying 
dynamics of the interaction between dark matter and nuclei. We will also show how 
anisotropies in the dark matter velocity distribution are manifest in rate anisotropies. 

A. Directional dependence from form factors and kinematics 

This subsection reviews the basic kinematics of dark matter direct detection. Consider 
a dark matter particle with mass m^m and initial velocity v^m scattering off of a nucleus at 
rest with mass m^- The angle between the nuclear recoil direction and Vdm, Vi * s given by: 

COST] = . (1) 

Vdm 

Above, v m i n (Efj) is the minimum dark matter velocity required by energy and momentum 
conservation in order for the nucleus to be scattered with recoil energy Er. If the kinematics 
of the scattering is elastic, v m m (Er) is given by 



= (2) 

where /i is the reduced mass of the dark matter-nucleus system. Fig. [T] illustrates the relation 
in Eq. [I] for different values of the initial dark matter velocity. 
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FIG. 1: Left: Kinematics of the DM- nucleus scattering. Right: The curves illustrate the 
relationship in Eq. [I]for different values of the initial DM velocity t>d m , assuming m dm = 100 GeV 
and the scattering is elastic off of a Germanium target. 

Notice that given v dra , the angle rj and the nuclear recoil energy E R are not uniquely 
determined. In fact, given v<im, the recoil energy can range in the interval 



< E R < 



m N 



(3) 



for elastic scattering. Given a flux of dark matter particles over the detector, the nuclear 
recoil rate is determined by Eq. [T] convolved with the dark matter velocity distribution and 
the scattering probability density over the energy recoil range in Eq. |3j If this probability 
distribution function (pdf) is flat in E R , the corresponding probability distribution over rj 
will range from to 7r/2, peaking at tt/4, as illustrated in Fig. [5J 

However, generically the scattering cross section will have a dependence on energy, 
modifying the probability that, given v^m, the nucleus will recoil with energy E R . For 
instance, the finite size of the nucleus will introduce a nuclear form factor dependence that 
suppresses the scattering probability at high recoil energies [37]. There may also be a form 
factor associated with the nuclear spin dependence of the interaction [38H4T)] , or on the dark 
matter properties, for instance if the dark matter has a dipole moment or is a composite 
state [HHIS]- Such form factors will alter not only the nuclear recoil energy distribution, 
but the nuclear recoil direction as well. 
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FIG. 2: Recoil energy (left) and recoil angle (right) distributions for a fixed dark matter velocity 
^dm- A Hat scattering probability over the kinematically allowed energy range will translate into a 
nuclear recoil directional pdf peaking at an angle of ir/4 with respect to the incident dark matter 
direction. Form factors will introduce an energy dependence of the scattering probability, distorting 
not only the recoil spectrum but also the nuclear directional rate. Kinematic properties such as 
inelasticity will not distort the pdf's, but truncate them over the kinematically allowed range. 

Fig. [2] illustrates how form factors distort the distributions on E R and rj for a fixed 
t>dm- For instance, the spin-independent nuclear form factor that parameterizes the loss of 
coherence in the scattering with the nucleus suppresses the probability of scattering at large 
recoil energies. The effect of that on directionality is to suppress the rate towards small rj. 
A dark matter form factor depending on a positive power of the momentum transfer, such 
as q 2 , will have the opposite effect, suppressing the scattering probability at low energies 
and shifting the recoil direction rate towards small angles with respect to the dark matter 
incident direction. 

Another interesting possibility is that the scattering is inelastic, and the dark matter 
particle up-scatters to an excited state with mass splitting 5 relative to the ground state 
In this case v m m(ER) will be given by 

I / m N E R 



Vmm(ER) 



y / 2m N E R V hn 



6). 



(4) 
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The presence of inelasticity will not alter the pdf 's shape in E R and rj for a given v<im, but 
only truncate the pdf to a narrower range [E Rmin , E RmeLX ] (and corresponding [r} m i n , r} max \) 
that satisfies the relation v dm = v min (E R )\ ERininjERm ^. 



B. Directional dependence from the velocity distribution 

As mentioned in the previous subsection, the energy dependence of the dark matter - 
nucleus interaction affects the directional rate of nuclear recoils. In this section we discuss 
the additional effect of anisotropics in the velocity distribution. 

So far we considered the effect of dark matter incidence with fixed velocity and direction. 
That scenario may be realized if the DM distribution in the solar neighborhood is dominated 
by a stream with low velocity dispersion (5UHS2]- Another possibility is that the dominant 
dark matter component in the solar neighborhood is isotropic and at equilibrium. As a first 
approximation we can consider a thermal velocity distribution truncated at some escape 
velocity, v e8c . If the dark matter halo and the solar system were co- rotating, this distribution 
would yield an isotropic nuclear recoil rate. However, if the solar system is moving through 
a static halo, the flux of dark matter particles will be boosted in the direction of the Earth's 
motion, creating interesting effects such as an annual modulation in the rate and a large 
anisotropy in the nuclear recoil direction, which will peak in the opposite direction to the 
Cygnus constellation. Anisotropic velocity distributions may generate even more interesting 
features in the directional rate of nuclear recoils, such as rings [53] and hot spots. 

In order to illustrate those features, we begin by establishing formulas and conventions 
that will be used throughout this paper. Given a dark matter velocity distribution f(v) in 
the galactic frame, the nuclear recoil rate as a function of energy and direction is given by 



dR 



N F{E R ) J d 3 v f(v) 6 [(v-v e ) ■ v R - v min (E R )\ (5) 



dE R dVt R 

where is the Earth's velocity in the galactic frame. Af contains the dependence on the 
scattering cross-section at zero momentum-transfer, on the local dark matter density and 
on phase-space factors. F(E R ) parameterizes the cross-section at finite momentum transfer 
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FIG. 3: (Left) The coordinate system adopted in this study. The x-direction points away from the 
center of the halo, and the z-direction is aligned with the axis of the disk. (Right) Angles 6 and cj) 
parameterizing the nuclear recoil direction: 9 is the angle between the vr and — u®, and 4> is the 
azimuthal angle with respect to the Earth's motion. 

as a product of nuclear and dark matter form factors. Finally, Vr is the direction of the 
nuclear recoil. We will ignore the Earth's orbit around the Sun and choose coordinates such 
that 

v® = ( , % , ), (6) 
v r = ( sin 9 cos cp , cos 9 , sin ^ sin ). (7) 

In those coordinates the ^-direction is aligned with the axis of the disk and the x-direction 
points away from the center of the galaxy as depicted in Fig. |3j We will also neglect the 
eccentricity of the Sun's orbit around the center of the galaxy. In that limit the Earth's 
velocity is parallel to the ^/-direction. 

Note from Eq. [5] that, if the velocity distribution f{v) is isotropic in the galactic frame, 
then the only source of anisotropy in the rate will come from the Earth's boost, £T e • -Or = 
% cos 9. Therefore, in the coordinate system we adopted, isotropic velocity distributions 
should induce rates that have a dependence on 9 but not on <ft. That means that any <p 
dependence on the rate is a smoking gun signature of an anisotropic velocity distribution. 
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FIG. 4: Nuclear recoil directions for dark matter particles scattering at high velocities (up) and 
low velocities (down) for isotropic (left) and anisotropic (right) dark matter velocity distributions 
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FIG. 5: 9 — <j> profiles of the directional detection rates for an isotropic dark matter velocity 
distribution at low recoil energy (left) and an anisotropic dark matter velocity distribution at low 
(middle) and high (right) recoil energies, where the anisotropy originates from suppressed circular 
orbits. 

For illustration purposes, we compare the angular dependence of the nuclear recoil rate 
for an isotropic versus an anisotropic distribution. As pointed out, the predicted rate for 
the isotropic distribution is flat on <ft, and has a dependence on 9 that varies with the recoil 
energy. If the energy is low, the recoil direction will prefer high values of 6, given Eq. [T]and 
the fact that v m i n (En) decreases as En decreases. On the other hand, increasing the recoil 
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energy will shift the rate towards lower values of 9. That is illustrated pictorially in Fig. Hi- 
low recoil energies exhibit "cone-like" features in the recoil direction that point away from 
the incoming dark matter flux. Increasing the recoil energy has the effect of collimating the 
"recoil cone", until the ring-like structure finally shrinks to a hot spot. 

Now consider an anisotropic velocity distribution for which the tangential components of 
the velocities are suppressed relative to the radial component. In that case circular orbits 
are subdominant, and the preferred flux of dark matter particles is along the radial direction 
in the galactic frame. In the Earth's frame, however, the flux picks up a component along 
the Earth's motion. That generates two overlapping "recoil cones" and hot spots in the 
recoil rate. As in the isotropic case, increasing the recoil energy will shrink the recoil cones 
to two more pronounced hot spots that will be shifted in by 7r/2 relative to the low energy 
hot spots. That effect is illustrated in Fig|5j 

III. INTEGRALS OF MOTION AND JEANS THEOREM 

In the previous section we discussed how the particle and astrophysical properties of dark 
matter manifest themselves in directional nuclear recoil signatures. Since very little is known 
about those properties, educated guesses have to made in order to explore direct detection 
signals. However, in the event that a dark matter signal is discovered in direct detection, 
we would like to take the data at face value and interpret it free of theoretical prejudice 
whenever possible. 

In this section we will develop the tools for the main goal of this work: to infer the 
dark matter phase space distribution in the solar neighborhood from directional data. The 
first step towards this goal is to choose the correct variables on which to express the phase 
space distribution function. Inquiring how the dark matter halo was formed in the first 
place will give us some guidance. For instance, since the Milky Way halo has not undergone 
recent mergers, we can assume that the majority of dark matter particles has had enough 
time to relax to a state of equilibrium, in which case Jeans Theorem should hold to a good 
approximation [58J: 

Jeans theorem : Any steady-state solution of the collisionless Boltzmann equation 
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depends on the phase-space coordinates only through integrals of motion. Conversely, any 
function of the integrals of motion is a steady-state solution of the collisionless Boltzmann 
equation. 

The Jeans theorem is the inspiration for the parameterization of the dark matter DF 
that will be used in this work. If a measurement of the local dark matter distribution 
function were our sole concern, any parameterization in terms of functions of the phase 
space coordinates could be adopted, even if such functions were not the true integrals of 
motion. But in scenarios where the equilibrium conditions of Jeans Theorem are valid, 
expressing the DF in terms of the true integrals of motion provides information beyond the 
scales of the solar neighborhood. In principle, we could analytically extend a measurement 
of the phase distribution in the solar system to any other position in the halo, and therefore 
infer the shape of the global DF through 

f(r,v) = f(I u I 2 ,...,I n ) = f{h{r,v)J 2 {r^,..,I n {r,v)l (8) 

1 1 . . . I n being integrals of motion of the halo. 

Moreover, the strong version of Jeans Theorem allows us to assume that, for all practical 

purposes, the phase-space distribution will be a function of at most 3 integrals of motion. 

Quasi-static haloes, for which Jeans theorem applies, always have energy, S, as an integral 

of motion. Other possibilities include the components of the angular momentum, L, and 

possibly non-classical integrals which have no analytical expression. Interestingly, the 

directional recoil rate is easily integrated if the distribution function depends only on energy. 

Consider a generic DF f(£). According to Eqj5j the directional rate is given by 

dR 
dE R dtt R 

where we have defined 



oc J 5(v-v R -p) d 3 v f(S) e(-£), (9) 



p = p(E R ,0,(j)) = v & -VR + v in i I1 (E R ), (10) 

and the energy 8 is given in terms of the dark matter's velocity and the local galactic escape 
velocity t> esc : 

£ = \-1>{?») = v —^. (11) 
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The factor 0(— S) excludes the possibility that the there are unbound particles (for which 
S > 0). 

On Appendix A we show that, 

5 (v-v R - P ) d s v G(-E) S" = 2n [^^j 6 (v 2 esc - p 2 ) . (12) 

Hence, expanding f(S) in a Taylor series, using Eq. 12 and re-summing the expansion we 
finally obtain: 

oc 2nQ{v 2 sc -p 2 ) f 2 f(S)dS. (14) 

J I v csc P \ 



Eq. 13 is a remarkable relation between the directional nuclear recoil rate and the 
underlying dark matter distribution. It makes it clear that, if the DF is isotropic in the solar 
neighborhood, measuring the dependence of the rate on v 2 sc /2 — p(E R , 9, (f>) 2 /2 is equivalent 
to directly measuring the underlying DF, since the later is the derivative of the rate with 
respect to v 2 sc /2 - p(E R , 9, (p) 2 /2. 

More generically, however, the dark matter DF will be anisotropic and depend on 
other integrals of motion, for instance the magnitude of the angular momentum \L\ or 
its components, such as L z . Unfortunately, no straightforward expression for the rate has 
been derived for a generic function f(S, \L\, L z ). However, we have been able to analytically 
integrate the rate for a certain class of functions (see Appendix [A]). Such functions have 
polynomial dependence in £ and harmonic dependence in the components of the angular 
momentum. That will prove useful for the method we will propose in the next section, where 
we will perform a special function decomposition of the DF. 



IV. SPECIAL FUNCTION DECOMPOSITION OF THE HALO DISTRIBUTION 

Although the reconstruction of the dark matter distribution function from directional 
and non-directional detection signals has already been investigated for specific halo types 
[32H3S], no model independent reconstruction study has been performed yet (see, however, 

13 



halo independent proposals to interpret direct detection signals [59H6T] ). One of the main 
goals of this paper is to provide a parameterization of the DF that is as generic as possible. 
The method we propose is to decompose the distribution function in an orthonormal basis of 
special functions, and fit the measured rate to the coefficients of this decomposition. If the 
basis is appropriately chosen, the higher order coefficients will be subdominant and we will 
be able to truncate the series to a finite and preferably manageable number of coefficients. 

There is no optimal choice for the variables and for the basis of decomposition that 
applies for all possible cases, and our main purpose is to illustrate out method rather than 
explore all possible variables and bases. Therefore, in this paper we will choose the following 
integrals of motion on which to parameterize the local dark matter distribution function: 

2 2 

• The energy S, which, as before, is given by £ = - ~ 2 " csc . 

• The ^-component of the angular momentum, L z , which in the coordinate system we 
have adopted is given by L z = \Jr 2 — z 2 v ( f > . In the solar neighborhood z = and hence 
L z = r^v^. 

• The magnitude of the angular momentum, whose DF dependence will be incorporated 
through the variable L t = ±\J \L\ 2 — L 2 Z = ^r 2 vj + z 2 v 2 . Again, in the solar 
neighborhood that reduces to L t = r^vg. 

We will make a further assumption which is not generic but that will greatly facilitate 
the computation of the rate and reduce the number of decomposition coefficients in the fit. 
This assumption is that the DF f(S, L t , L z ) is separable: 

f(£,L t ,L z ) = f x (E)f 2 (L t )h{L z ). (15) 

Under this assumption, we will adopt the following choice of basis for fi(£), f2i.Lt) and 
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ML,): 

f 2 {L t ) = J^c^cosfnn-^-), (16) 

n \ L max / 

h{Lz) = ^4 m cos (mvr-^- J . 

Above, is expanded in shifted Legendre polinomials, Pi, and £i im is the lowest dark 

matter energy that an experiment with recoil threshold is able to probe: 

_ K in (£g) - %) 2 - ^ sc ^ 

In principle we could have chosen £ lim to the lowest possible energy, —v1 sc /2. However, since 
a realistic experiment will have no information about particles with energy less than £\i m , 
there is no purpose in parameterizing the distribution function beyond what is measurable. 
The functions f^Lt) and fs{L y ) are decomposed in a Fourier series, where L max = r e f esc - 

A great practical advantage of the parameterization above is that the recoil rate can 
be analytically integrated. That way no numerical integration is necessary, reducing 
computational time to perform the fits. In Appendix [A] we derive and list the analytical 
forms for the rate on those basis functions. 

In the following section we will illustrate how the fit is performed for a hypothetical 
experiment measuring the recoil rate for an anisotropic DF. We will choose as an example 
a Michie distribution [U2] and generate mock up Monte Carlo data as the detector would 



measure it. We will then employ our method to recover the original distribution. In Sec. VI 



we will use the Via Lactea II simulation as the underlying dark matter distribution and 
simulate the rate at the detector. We will then employ our method to perform the fits. 



V. DIRECTIONAL FITS TO A MICHIE DISTRIBUTION 



In this section we illustrate our method for a hypothetical dark matter signal in a 
directional detector. Ideally we would like to have a large number of events to perform the 
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fit. However, the current limits are quite stringent for spin-independent elastic scattering 
dark matter particles with mass of (9(100) GeV. For instance, the latest XENON100 run 
with 4 kg-yrs of exposure observed 3 candidate events with an expected background of 
1.8 ± 0.6 events [21] . That suggests that, in the heavy WIMP scenario where discovery is 
just around the corner, the detection of (9(100) events would require at least hundreds of 
kg-yrs of exposure, which is still an unrealistic goal for the (D(l)-kg detectors planned by 
current directional experiments. The most optimistic possibility is if dark matter is light, 
of order of a few GeV. In this mass range the limits are much less stringent [H31 El], and 
the relatively light targets considered by directional experiments, such as CF 4 and CS 2 , are 
ideal. 

Therefore we choose to generate this hypothetical signal for an elastically scattering dark 
matter particle with mass md m = 6 GeV. We assume that the scattering is off of a CS2 
target, with recoil energy acceptance window Er G [5keV, 15keV] and for simplicity we 
ignore the scattering off of carbon. We assume that the experiment operates in a zero 
background environment, but in principle background events could be straightforwardly 
incorporated in the fits if their directional shape is known. Finally, the results we will show 
in the following are obtained from fits without resolution effects taken into account. As a 
crude approximation of the detector angular resolution, we have divided the nuclear recoil 
direction in bins of size (10°, 10°) and did not observe any significant alteration of our results. 
In the results shown in the rest of this article, no binning has been used. 

In order to illustrate the interesting effects of anisotropy, we choose a Michie distribution 
for the dark matter. In terms of integrals of motion, the Michie distribution is given by: 

/Michie^, \L\)= N e -^l 2 / L ° (e- 2£ ^-l) Q{-£) (18) 

Our choice of parameters is the following: vq = 280 km/s, v esc = 600 km/s, L = r©t> and 
a = 1. We will assume in the fits that local escape velocity f esc as well as the dark matter 
mass mdm are know to sufficient precision. That is a reasonable assumption since directional 
experiments will only come into play after the recoil spectrum has been well measured by 
several non-directional experiments with different nuclear targets. 



Note that the Michie distribution (18) trivially satisfies the separability condition (15) 
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when expressed in terms of S, L t and L z , since \L\ 2 = L\ + L z . Then, with 

fi(£) = e-^-l, (19) 
f 2 (L t ) = (20) 
/ 8 (L.) = e- QL ^ L o (21) 

we can perform the special function decomposition as in EqJTBJand obtain the coefficients: 

{c P0 , c Pl , c P2 , c F3 , C P4 ,CP5,...} = {1,1.65,0.890,0.297,0.0722,0.0137,...}, (22) 
{c F0 ,c F1 ,c F2 ,c F3 ,...} = {0.412,0.485,0.0952,0.00773,...}, (23) 
{c F0 ,c z F1 ,c z F2 ,c z F3 ,...} = {0.412,0.485,0.0952,0.00773,...}. (24) 

The overall normalization of the coefficients above is unphysical since it can be reabsorbed in 
the cross section. Hence, the coefficients were scaled such that fi(£), f^{L t ) and fz{L z ) are 
normalized to unity. We have truncated the energy expansion in Legendre polynomials to 
the 5th order, and the Fourier decomposition of the components of the angular momentum 
to the 3rd harmonic. The truncation of the expansions at this orders approximates the true 
functions to an accuracy better than the percent level. 

Knowing the true values of the decomposition coefficients, we can compare them to the 
values obtained from fits of the rate. We generated several Monte Carlo data ensembles, 
varying the number of observed events. We then performed unbinned log-likelihood fits in 
the 11-dimensional parameter space of decomposition coefficients using the publicly available 
code MultiNest v2.10, a Bayesian inference tool that uses a nested sampling algorithm to 
scan the parameter space [651 EE] • 

The 1-cr range of the fitted coefficients can be extracted by approximating the log- 
likelihood function near its maximum by a quadratic function of the parameters of the 
fit: 

C{C)^C m3X - l -Y.°l^ ( 25 ) 

i 1 

—* 

where C is the vector representing a given set of coefficients and a is their standard deviation. 
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FIG. 6: Results of the fits for fi(£), /bC^t) an d fc(L z ), assuming 10 4 (red) and 10 3 (blue) signal 
events. The bands delimit the ranges of these functions given by the 1-er errors on the decomposition 



coefficients. The underlying phase space distribution is given by the Michie DF in (18), represented 
in graphs above by the solid lines. 

In Fig. [7] we show the results of the fit for different numbers of events, varying from 100 to 
10,000. We display the 68% C.L. contour intervals of the fitted coefficients, as well as the true 
values for comparison. As the errors are correlated, the 68% range of the decomposition 
coefficients is actually a multidimensional ellipsoid. The different plots in Fig. [7] are 2D 
projections of this ellipsoid in different planes of the parameter space. 

As one can see, a few thousand events are necessary for a good measurement of the 
underlying dark matter distribution. Fig. [6] displays the functional forms of fi(S), f%{L t ) and 
fz{L z ) resulting from fits with 1000 events and 10,000 events. This shows that our method 
works remarkably well if the variables and decomposition basis are chosen judiciously. In 
the real world, however, an appropriate choice of decomposition basis may not be as trivial 
as in this hypothetical example. However, several checks of the data can be made before 
performing fits to test whether a given hypothesis is a good approximation. We will illustrate 
that in the next section, where we will use the Via Lactea II N-body simulation as a template 
for the dark matter distribution in our galaxy, and "scatter" the particles of this simulation 
is a detector placed at a location around 8 kpc from the center to try to infer the underlying 
phase space distribution. 
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FIG. 7: 68% confidence contours on the decomposition coefficients of the DM distribution, found by 
fitting several mock directional rates generated by Monte Carlo, with different numbers of observed 
events. The dots correspond to the true values of the coefficients. 

VI. FITS TO THE VIA LACTEA II SIMULATION 



No dark matter direct detection signal has been observed to this date to provide us 
with clues on the local distribution. Most of our current understanding of the range of 



19 



possibilities for the local distribution comes from N-body simulations [U |50j EH [67], which 
also provide one of the best test grounds to investigate the formation and dynamics of dark 
matter halos [BE]- With about 10 9 particles with masses around 10 3 solar masses, Via Lactea 
II is currently one of the highest resolution N-body simulations and can be used as a model 
of the Milk Way halo [69]. 

In this section we will test the method proposed in this study using Via Lactea II as 
the underlying dark matter distribution. Ideally, we would like to choose a box located at 
= 8 kpc away from the center of the halo (which is the distance of the solar system 
from the center of the Milky Way) and with a volume (Ar) 3 <C where Ar is the length 
of a side of the box. Unfortunately, in Via Lactea II a region of this size typically does 
not enclose any particles. Our approach is then to choose a spherical shell of 8kpc-radius 
and 200pc-thickness, which encloses enough particles to generate a dark matter signal at 
a directional detector. Since we do not want to wash out anisotropies, we choose a fixed 
position r e at which to place the detector, and for each dark matter particle in the shell 
with coordinates v and r, we perform an affine map v{r) — > v(r^). 

In the same spirit of the previous section, we would like to fit a signal generated from Via 



Lactea II to a distribution function in terms of £ , L t , L z , which is separable as in Eq. [T5j and 
decomposed as a Legendre series in £ and a Fourier series in L t , L z , see Eq. [16} However, 
before attempting to perform those fits we should test the validity of those assumptions. 

The first of those assumptions, that the energy and angular momentum are integrals of 
motion, can be shown to be a good approximation for positions not too far away from the 
center of the halo, where the gravitational potential can be reasonably approximated as 
being spherical. A quantitative measure of the sphericity of the potential is given by 

i>{f)- ip s (r) 



= ^ > ™ > (26) 

Above, ip(f) is the exact gravitational potential at the position r, and ip s (r) is the 
potential computed assuming spherical symmetry and evaluated using Gauss' law. x(r) 
then parametrizes the deviation of the gravitational potential at position r*from a spherically 
symmetric potential. 

Fig. [8] displays the r dependence of \ h 1 three orthogonal directions: r = rx, f = ry 
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FIG. 8: The deviation of the gravitational potential from spherical symmetry, parameterized by 



x(r), see Eq. 26 Three orthogonal directions for r are displayed: r = rx (blue), f = ry (green) 
and f = rz (red). 

and f = rz. One can see that for r < 200 kpc, the assumption of a spherically symmetric 
potential leads to an error (9(6%). For r < 30 kpc, the error is reduced to 0(3%). Hence 
the energy and angular momentum of dark matter particles within a 30 kpc radius should 
be approximately conserved. 

The second of our assumptions, that of separability, can be evaluated through the function 

g 2 f(£,L t ,L z ) 



G(£, L t , L z ) 



g e (£)g Lt {L t )g Lz {L z y 



where 



9e(S) = j f(S,L t ,L z )dL t dL 2 



9L t {L t ) 



f(£,L t ,L z )d£dL z 
f(£,L t ,L z )d£dL t 



g = J f(£,L t ,L z )d£dL t dL z 



(27) 

(28) 
(29) 
(30) 
(31) 



which measures the correlation between the three variables £, L t and L z . If the dark matter 



distribution function verifies Eq. 15, G(£,L t ,L z ) should be constant and equal to 1. 
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FIG. 9: Gg, Gt t and Gl. in function of £, Lt and L z , respectively. 



Due to the lack of statistics, in Fig. |9j we plot Gg, Gl v , Gl z , the averages of G over 
the two other variables, rather than the tridimensional behavior of the function. The fact 
that Gg, GrL t , Gl z are consistent with a constant value equal to one reassures us that our 
assumption of separability is well justified. 



A. Local Fits to Via Lactea II 



The range of validity of our hypotheses described in Sec. [ITTjwere quantified in the previous 
section. In this section we perform fits to a directional detection signal generated from the 
distribution of VLII simulation, using particles contained within a shell located 8 kpc away 
from the center of the simulation, and with thickness of 200 pc. As mentioned previously, 
we fix the location for the detector at r©, and transform the velocities of all particles in the 
shell through an affine map v(f) — > v (?"©). 

As in Sec [V] we set the mass of the dark matter particles to = 6 GeV. Inspection 
of the VLII distributions allows us determine that the escape velocity at 8 kpc is f esc ~ 
545 km/s. We generate 10,000 events by elastically scattering the particles in the shell off 
of a CS2 target, but ignore the scattering off of carbon. However, differently from Sec [Vj 
we are forced to lower the energy threshold of 5 keV. The reason for that is that with a 
5 keV threshold, only particles in the tail of the energy distribution (v > 518 km/s) are 
kinematically allowed to scatter. Unfortunately there only about 10 5 particles inside the 
shell, and the tails of the distributions of those particles are not well populated. Hence, in 
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order to avoid large errors resulting from fitting the tails, we lower the energy threshold to 
2 keV, allowing particles with v > 330 km/s to scatter. In that range the distribution is 
better populated and the fits are more reliable. As in Sec [Vj detector resolution effects are 
not taken into account. 
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FIG. 10: 6 — (j) profile of the directional signal generated from VLII simulation. Notice that the 
hint of hot spots at (f> = [tt] suggests the presence of anisotropics. 



The 9 — (p profile of the directional detection rate generated is shown in Fig. [TOj The hint 
of hot spots at = [tt] suggest a suppression of circular orbits perpendicular to the plane 
of the Earth's orbit in the halo. At this point this is just a qualitative statement, which will 
be made quantitative below by the results of the fits. Note that the plane of the Earth's 
orbit was chosen arbitrarily given that VLII does not have a baryonic disk. 

We perform unbinned log-likelihood fits to this directional signal as described in 
Sees. IVjV The results of the fits are displayed in Fig. 11, where we plot £, L t and L z 



distributions for the particles that are kinematically allowed to generate nuclear recoils 
within the recoil energy acceptance window. We contrast the distributions extracted from 
the fit with the true underlying distributions to show that the fit is correctly capturing the 
local VLII phase space distribution within errors. The fits for the £, L t and L z dependent 
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FIG. 11: Results of fits to a directional signal generated with VLII, displayed as distributions of 
£ , L y and L z for the particles kinematically allowed to scatter. The red band corresponds to the 
fitted distributions and its 1 a error, whereas the blue histograms correspond the true underlying 
distributions. 



parts of the DF, denoted by fi(£), /^(-^t) and fz(L z ) in eq. 15, are displayed as blue bands 



in Fig. 13. One can see that although fz(L z ) is consistent with a flat function, f2(L t ) is 
suppressed for high values of L t , meaning that orbits perpendicular to the plane of the 
Earth's orbit in the halo are disfavored. That confirms the presence of anisotropy mentioned 
earlier, and is consistent with the hot spots developed in the rate at <fi = [tt], see Fig. 10 



B. Global Extrapolation of the Local Fits 



As discussed in Sec. |III[ Jeans Theorem tells us if the dark matter distribution function 
can be expressed locally in terms of integrals of motion, then it should be valid at all other 
positions of the halo, provided that the distribution is in a steady state of equilibrium. In 
this section we test the validity of our choice of £, L t and L z as integrals of motion by 
checking whether the VLII distribution function at other positions in the simulation agrees 
with our fit at 8 kpc. 

The most straightforward variable to check is £. Consider the distribution of dark matter 
particles at a given location r, binned in £, L t and L z with bin widths A£, AL t and AL Z 
respectively. The number of particles in the bin centered around £, L t and L z is given by: 

N(S, L t , L z ) = A£AL t AL z Det[J(£, L t , L z )\ f x (£)f 2 {L t )f 3 (L z ), (32) 
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FIG. 12: fi(£) extracted from VLII data for r ~ 8 kpc (blue) and r ~ 30 kpc (green), using Eq. 35 
The values of £ are normalized to £u m (r) = —ip(r) at r = 8 kpc. 

where Det[J(£, L y , L z )] is the determinant of the Jacobian, given by: 

Det[J(£, L t , L,)] = 1 2 ^tT ( 33 ) 

V 2S + ' u csc - L t/ r2 - L2 z/ r2 

Since the bins with better statistics are those with lower values of £, and hence lower 

values of L t and L z , we can extract f\(£) from bins where L t and L z are centered around 

zero and £ is centered around £f. 

f(£)= 1 1 ^.o.o) (M) 

Jl[t) / 2 (0)/ 8 (0) A£AL t AL z J 2 £, + vl r ' 1 J 



The first factors in Eq. 34 are independent of £ and can be normalized away, finally 
giving: 



fi{£i) 



N(£ h 0,0) 



(35) 



In Fig. 12, we plot fi(£) using Eq. 35 for two different positions in the VLII simulation, 
namely r = 8 kpc and r = 30 kpc. The distribution of particles at different radii are 
obtained from shells with the corresponding radii and 200 pc-thickness. Note that since the 
gravitational potential decreases with the distance to the center of the galaxy, the minimum 
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FIG. 13: Fits for (a) (b) f 2 (L t ), and (c) f 3 (L z ) at 4.5 kpc (red), 8 kpc (blue) and 30 kpc 

(green) using the distribution of Via Lactea II. The black dots with error bars in (a) correspond to 



values of fi{£) extracted directly from data at 8 kpc using Eq. 35 The good agreement at different 
radii corroborates the use of Jeans Theorem. 

energy of the particles at r = 30 kpc (£ m m( r ) = ~ ^( r )) is higher than at r = 8 kpc. The 
excellent agreement between the extraction of fi{£) for the two radii reinforces the validity 
of Jeans Theorem and the fact the £ is a good integral of motion. 

Unfortunately, there is not enough statistics to extract f2i.Lt) and fz{L z ) directly from 
the simulation using the same method used above for fi(£), since the population in bins 
with increasing values of L t and L z rapidly drop in number. The alternative is to fit the 



dark matter distributions at different radii using our parameterization in Eqs. 15 and 16 



We perform such fits for three different radii, r = 4.5 kpc, r = 8 kpc and r = 30 kpc, where 
as before we select particles within shells of 200 pc-thickness and radius r. 



Fig. 12 shows the results of the fits for fi(£), f2(Lt) and f 3 {L z ) and their corresponding 
error bands. The fact that there is a reasonable agreement over a wide range range of 
galactic radii is evidence that the parameterization motivated by Jeans theorem is valid and 
meaningful. From Fig. [12] it is also interesting to note that, although fz{L z ) is consistent 
with a flat distribution, /^(-^t) indicates that the VLII distribution is anisotropic, and the 
angular momentum component along the Earth's motion is suppressed. Another interesting 
fact is that f\{£) is not consistent with a single exponential, as is the case for Maxwellian 
or Michie distributions. Indeed, the high velocity tail of fi{£) is harder than Maxwellian, 
as pointed out in [57] . 
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VII. DISCUSSION 



In this paper we discussed important aspects of directional detection that will be relevant 
once dark matter is detected in underground experiments. Directional experiments will 
probe the dynamics and kinematics in the dark matter interactions with nuclei and allow us 
for directly measure the dark matter phase space distribution in the solar neighborhood. 

It will be important to reduce our bias in extracting information from directional data. 
With that in mind, we proposed to parametrize the local distribution function in terms of 
integrals of motion. If the dominant component of underlying dark matter distribution in 
the solar neighborhood is in equilibrium, expressing the DF in terms of integrals of motion 
will allow us to infer the DF at other positions in the halo. That can have a potentially large 
impact in studies of indirect signals, as well as implications for our understanding of the 
formation of the Milky Way halo. There is a chance that the local dark matter distribution 
could have a sizable component that is in a non-equilibrium configuration such as streams 
of dark matter arising from tidally disrupted dark matter sub-halos. These features can be 
readily identified in directional detection experiments and can be subtracted off [5UH52]- 

We also proposed a decomposition of the distribution function in special functions of 
integrals of motion. That allows for a more model independent extraction of the true 
phase space distribution. We illustrated our method by generating Monte Carlo signals 
in a hypothetical directional detector and fitting the coefficients of our decomposition. This 
study was performed with an analytical model of the halo, as well as using Via-Lactea II 
simulation. In the later, we also tested the results of our fit by extrapolating the fitted local 
distribution to other positions in the halo and comparing to the true distribution. We found 
that for radii less than ~ 30 kpc our extrapolation worked to a very good approximation. 

The number of signal events necessary to detect non-standard features of the dark matter 
distribution obviously depend on the degree of such effects. For Via Lactea II we found 
that 0(1000) events were necessary to detect presence anisotropics and departures from a 
Maxwellian distribution. 

The parameterization and decomposition of the dark matter phase space distribution we 
have proposed may also have applications beyond directional detection. For instance, it could 
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also allow a better understanding of the dark matter haloes obtained in N-body simulations. 
In particular, reconstructing the dark matter phase space distribution would allow us to 
compare the different existing N-body simulations to each other, studying for example how 
the anisotropies of the dark matter distribution function depend on the different parameters 
of the simulation, especially the presence of a baryonic disk. 

In conclusion, this article has shown the power of directional detection experiments in the 
post-dark matter-discovery era. By using the results from directional detection experiments 
and N-body simulations, a better understanding of dark matter in the Milky Way can be 
obtained. 
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Appendix A: Analytic Scattering Rates 

In this appendix we derive useful analytical expressions for the directional rate of nuclear 
recoils. As discussed in Sec. [TT], the differential nuclear recoil rate is given by 
dR 



dE R dQ 

where 



oc 



J dvf(v) 5[v-v R -p(E R ,e)]G[v 2 esc -v 2 }, (Al) 



p(E R , 9) = v min (E R ) -v S) -v R = v min (E R ) - v e cos 9. (A2) 
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Using 



v r — ( x r, Vr-i z r) = (sin 9 cos 0, cos 9, sin 9 sin 0), 



and integrating over v x gives 
dR 



oc 



dE R dVt \x R 



^— [ dv y dv z f(v x [v y , v z ],v y , v z )Q[v 2 sc - v x [v y , v z } 2 -v\- v 2 z ], (A3) 

>R\ J 



where 



v x (v y ,v z ) = —(p(E R ,6) - y R v y - z R v z ). 

Xr 

The domain of integration for the variables v y and v z is then set by the condition 



(A4) 



v z (zr + x r) + 2v z z R (y R v v - p(E R , 9)) + (p(E R , 9) - y R v y ) 2 - x 2 R (v 2 sc - v 2 ) < 0. (A5) 



The limits of integration on v z are then 



Vz±(Vy) = 



X R + Z R 



z R {p{E R , 9) - y R v y ) ± VA 



where 



a = 4 [(4 + 4)( 



2 2 



)-(p(E R: 9)-v y y R ) 2 ]. 



The condition A > sets the limits for v y : 



v y ± = y RP (E R , 9)±y/(l- y R )(vl c - p\E R , 9)). 
With the limits of integration defined, the only remaining constraint is 

v 2 sc -p 2 (E R ,9))>0, 



(A6) 



(AT) 



(A8) 



(A9) 



which requires that for a given recoil energy E R and recoil angle 9, the minimum velocity 
of the incoming dark matter particles in the galactic frame must be smaller than the dark 
matter escape velocity. 

We then do the variable change 



V y = v y - y R p{E_ 



R, 



Vy ~ V y , 



v z = X 1±A, 

\xr\ 



: ^[p(E R ,9)(l-y 2 R )-y R Vy] = Xl1 [ 1 



zr 



-v, - v". 



xr 



xr 



(A10) 
(AH) 
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The new boundaries are then given by ±V®(V y ) and ±V® with 

V z °(V y ) = ^Al - VI 



where 



V y = A o, 



AZ = (l-V&(vL-ftE R ,0)). 



(A12) 
(A13) 

(A14) 



The differential rate then becomes 



dR 



oc 



dE R dn x R + z 2 R J_ Ao J_ v ^Zv2 



A ry/^l-V* 



Q[v 2 esc -p 2 (E R ,9)} dV z dV y 



v x [V y + vl 



\ X R\ 



y ^ u y „2 



X R + Z R 



(y z + v z )),V y + v° y 



\ x r\ 



y „2 



X R + Z R 



(A15) 



For a distribution function expressed in terms of energy £ and the components of the 

angular momentum L t and L z , we can recover the dependence on V y and V z using 

1 
2 



£ = =;(v y 2 + v 2 -Al) 

L t = 

L z = r^Vy + v. 



J-' ,v> I f^rp 



(A16) 

(A17) 
(A18) 



Eqs. A15| to A18 allow us to integrate the rate analytically for several special forms of 
f(£,L t ,L z ), such as 



£ cos(mrL z ) cos(mTcL t ), 
E l cos(mrL z ) sm(mirL t ), 
£ e sm(nnL z ) cos(m-RL t ), 
£ e sm(nnL z ) sm(mnL t ). 



(A19) 
(A20) 
(A21) 
(A22) 



For instance, for /(£, L t , L z ) as in Eq. A19, the differential rate is given by 
dR 1 



1 /-Ao r^=^l 

dE R dn ^2'(x| + 4) y_ Ao L^=vi 



oc- 



cos 



:(V z + v° z ) 



(A23) 



cos 



[mrr 9 (V y + v°)] dV z dV y . 
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The integration above is straightforward to perform and we finally obtain: 
J 5 (v ■ vr — p) d 3 v E l cos(nirL z ) cos(mivL t ) = 



{-if t\ 7r cos[np(ny R + mz R )] x M + 
+(— 1)* ft 7T cos[irp(ny R - mz R )\ x M_, 



where we set r e = 1 to lighten the notation and 
p = v E cos 9 + v min (E R ), 

v R = (x R ,y R , z R ) = (sin 9 cos 0, cos 9, sin 9 sin 0), 



/ 2 _ 2 \ m 

m + = ( 9 , = , 2 %(W»» 2 + ™ 2 - Nto + mz R y^i^), 



^TTA/n 2 + m 2 - (ra/^ - mz R ) 2 J 
Similar expressions hold for the others: 

J 5 (v • v R — p) d 3 v E l sm(nirL z ) cos(mnL t ) = 

(-if l\ 7r sm[np(ny R + mz R )\ x M + 
+(-1)' l\ vr s\n[itp(ny R - mz R )\ x M_, 



J 8 (v • v R — p) d 3 v E l cos(n7vL z ) s'm(m7vL t ) = 

{-if l\ 7T sm[np(ny R + mz R )\ x M + 
+(— 1)* £\ ii sm[irp(-ny R + mz R )\ x M_, 
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J 5 (v • v R — p) d 3 v E l sin(n7rL z ) sin(m7rL t ) = 

(-1)' t\ n (- cos[7r p(ny R + mz R )}) x M H 
+(— l) e l\ it cos[7rp(ny R — mz R )\ x M_. 

Note that all those expressions can be written in a compact form: 
J 5 (v ■ v R - p) d 3 v £ e e iw{nL * +mLt) = (-l) e l\ 2n e ™p(nyR+mz R ) 

V V csc - P 



x . = ; Y +1 -WW" 2 + » 2 - (ny R + ^rYV^?) 

7iyn z + m z — [ny R + mz R ) z 



Remembering that 

vr = (x R: y R ,z R ), 
L = (0,L z ,L t ), 

we can rewrite the expression above in a frame independent form: 



J ' S(v-v R -p) d 3 v £ e e il ^ T = (-l) e i\ 2tt e^ T (v 2 csc - p 2 ) e+1 



with E = ^/v 2 sc - p 2 ^/l - (v R -u T ) 2 , 

where ut is an arbitrary vector in the tangential plane. 

Finally, when there is no dependence on L the rate takes the form: 



Je+i(\ur\ 




(\ut\ 
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In particular, for a generic /(£), 

J S(v-v R -p)d 3 v f(£) 
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